//
// GPSMapEdit
// (c) Konstantin Galichsky (kg@geopainting.com), 2002-2007
//
// Geodetic datum handling.
//

# include "StdAfx.h"
# include "Datum.h"
# include "Ellipsoid.h"
# include "PtAlgo.h"
# include "Globals.h"

const CEllipsoid c_eAiry1830        ("Airy 1830",               6377563.396,  1/299.3249646);
const CEllipsoid c_eAiryMod         ("Modified Airy",           6377340.189,  1/299.3249646);
const CEllipsoid c_eAustralian      ("Australian National",     6378160.0,    1/298.25);
const CEllipsoid c_eBessel1841      ("Bessel 1841",             6377397.155,  1/299.1528128);
const CEllipsoid c_eBesselNamibia   ("Bessel 1841 (Namibia)",   6377483.865,  1/299.1528128);
const CEllipsoid c_eBesselNorway    ("Bessel 1841 (Norway)",    6377492.0176, 1/299.1528);
const CEllipsoid c_eClarke1858      ("Clarke 1858",             6378350.87,   1/294.26);
const CEllipsoid c_eClarke1866      ("Clarke 1866",             6378206.4,    1/294.9786982);
const CEllipsoid c_eClarke1880      ("Clarke 1880",             6378249.145,  1/293.465);
const CEllipsoid c_eClarkeIGN       ("Clarke 1880 IGN",         6378249.2,    1/293.466021);
const CEllipsoid c_eClarkePalestine ("Clarke 1880 Palestine",   6378300.789,  1/293.466);
const CEllipsoid c_eEverest1948     ("Everest (1948)",          6377304.063,  1/300.8017);
const CEllipsoid c_eEverestIndia    ("Everest (India 1830)",    6377276.345,  1/300.8017);
const CEllipsoid c_eEverestIndia56  ("Everest (India 1956)",    6377301.243,  1/300.8017);
const CEllipsoid c_eEverestMalasia  ("Everest (Malaysia 1969)", 6377295.664,  1/300.8017);
const CEllipsoid c_eEverestPakistan ("Everest (Pakistan)",      6377309.613,  1/300.8017);
const CEllipsoid c_eEverestSS       ("Everest (Sabah Sarawak)", 6377298.556,  1/300.8017);
const CEllipsoid c_eEverestTimbalai ("Everest (Timbalai)",      6377298.561,  1/300.8017);
const CEllipsoid c_eFischer1960     ("Modified Fischer 1960",   6378155.0,    1/298.3);
const CEllipsoid c_eGRS80           ("GRS 80",                  6378137.0,    1/298.257222101);
const CEllipsoid c_eHayford         ("Hayford 1909",            6378388.0,    1/296.959263);
const CEllipsoid c_eHelmert         ("Helmert 1906",            6378200.0,    1/298.3);
const CEllipsoid c_eHough           ("Hough 1960",              6378270.0,    1/297.0);
const CEllipsoid c_eIndonesian      ("Indonesian 1974",         6378160.0,    1/298.247);
const CEllipsoid c_eInternational   ("International 1924",      6378388.0,    1/297.0);
const CEllipsoid c_eKrasovsky1940   ("Krasovsky 1940",          6378245.0,    1/298.3);
const CEllipsoid c_eSA1969          ("South American 1969",     6378160.0,    1/298.25);
const CEllipsoid c_eWGS72           ("WGS 72",                  6378135.0,    1/298.26);
const CEllipsoid c_eWGS84           ("WGS 84",                  6378137.0,    1/298.257223563);

///////////////////////////////////////////////////////

	// Arc seconds per radian.
const float c_ro = 206264.8062f;

CMolodenskyTransformation::CMolodenskyTransformation (
	const CEllipsoid & _SourceEllipsoid, const CEllipsoid & _DestEllipsoid,
	float _fDX, float _fDY, float _fDZ,
	float _fElevation
) :
	dx (_fDX), dy (_fDY), dz (_fDZ),
	fElevation (_fElevation),
	a  ((_SourceEllipsoid.a  + _DestEllipsoid.a)/2),
	e2 ((_SourceEllipsoid.e2 + _DestEllipsoid.e2)/2),
	da  (_DestEllipsoid.a  - _SourceEllipsoid.a),
	de2 (_DestEllipsoid.e2 - _SourceEllipsoid.e2)
{}

void CMolodenskyTransformation::ApplyTo (point_t & _point) const {
	const float & Ld = _point.x;
	const float & Bd = _point.y;

	const float B = Bd*c_PI/180;
	const float L = Ld*c_PI/180;

	const float fSinB  = ::sin (B);
	const float fCosB  = ::cos (B);
	const float fSinL  = ::sin (L);
	const float fCosL  = ::cos (L);

	const float fNu = 1 - e2*fSinB*fSinB;
	const float fSqrtNu = ::sqrt (fNu);

	const float M = a*(1 - e2)/(fNu*fSqrtNu);
	const float N = a/fSqrtNu;

	//
	// Latitude.
	//
	const float dB = c_ro/(M + fElevation)*(
		N/a*e2*fSinB*fCosB*da
		+ (N*N/(a*a) + 1)*N*fSinB*fCosB*de2/2
		- (dx*fCosL + dy*fSinL)*fSinB
		+ dz*fCosB
	);
	const float fWGS84Lat = Bd + dB/3600;

	//
	// Longitude.
	//
	const float dL = c_ro/((N + fElevation)*fCosB)*(-dx*fSinL + dy*fCosL);
	const float fWGS84Long = Ld + dL/3600;

	_point.x = fWGS84Long;
	_point.y = fWGS84Lat;
}

///////////////////////////////////////////////////////

CHelmertTransformation::CHelmertTransformation (
	const CEllipsoid & _SourceEllipsoid, const CEllipsoid & _DestEllipsoid,
	float _fDX, float _fDY, float _fDZ,
	float _fWX, float _fWY, float _fWZ,
	float _fMS,
	float _fElevation
) :
	CMolodenskyTransformation (_SourceEllipsoid, _DestEllipsoid, _fDX, _fDY, _fDZ, _fElevation),
	wx (_fWX), wy (_fWY), wz (_fWZ),
	ms (_fMS)
{}

void CHelmertTransformation::ApplyTo (point_t & _point) const {
	const float & Ld = _point.x;
	const float & Bd = _point.y;

	const float B = Bd*c_PI/180;
	const float L = Ld*c_PI/180;

	const float fSinB  = ::sin (B);
	const float fCosB  = ::cos (B);
	const float fTanB  = fSinB/fCosB;
	const float fCos2B = 2*fCosB*fCosB - 1;

	const float fSinL  = ::sin (L);
	const float fCosL  = ::cos (L);

	const float fNu = 1 - e2*fSinB*fSinB;
	const float fSqrtNu = ::sqrt (fNu);

	const float M = a*(1 - e2)/(fNu*fSqrtNu);
	const float N = a/fSqrtNu;

	//
	// Latitude.
	//
	const float dB = c_ro/(M + fElevation)*(
	    N/a*e2*fSinB*fCosB*da
		+ (N*N/(a*a) + 1)*N*fSinB*fCosB*de2/2
		- (dx*fCosL + dy*fSinL)*fSinB
		+ dz*fCosB
	)
		- wx*fSinL*(1 + e2*fCos2B)
		+ wy*fCosL*(1 + e2*fCos2B)
		- c_ro*ms*e2*fSinB*fCosB;

	const float fWGS84Lat = Bd + dB/3600;

	//
	// Longitude.
	//
	const float dL = c_ro/((N + fElevation)*fCosB)*(-dx*fSinL + dy*fCosL)
		+ fTanB*(1 - e2)*(wx*fCosL + wy*fSinL) - wz;

	const float fWGS84Long = Ld + dL/3600;
/*
	//
	// Elevation.
	//

	const float dH = -a/N*da + N*fSinB*fSinB*c_de2/2
		+ (dx*fCosL + dy*fSinL)*fCosB + dz*fSinB
		- N*e2*fSinB*fCosB*(wx/c_ro*fSinL - wy/c_ro*fCosL)
		+ (a*a/N + _fElevation)*ms;

	const float fWGS84Elevation = fElevation + dH;
*/
	_point.x = fWGS84Long;
	_point.y = fWGS84Lat;
}

///////////////////////////////////////////////

# ifdef _DEBUG
void Test_Datums () {
	point_t p (37.801902f, 55.6f);
	const CTransformationBase * p42 = NULL;
	FindOziDatumToWGS84 ("Pulkovo 1942 (1)", p42);
	p42->ApplyTo (p);
	assert (::fabs (p.x - 37.8f) < 5e-5);
	assert (::fabs (p.y - 55.6f) < 5e-5);
}
# endif // _DEBUG

///////////////////////////////////////////////
// New datums can be added using information on http://earth-info.nima.mil/GandG/tr.html
// "DoD World Geodetic System 1984, Its Definition and Relationships with Local Geodetic Systems"

CDatum g_Datums [] = {
	{"Adindan",                       "Adindan",                            & c_eClarke1880,      -162,  -12,  206},
	{"Afgooye",	                      "Afgooye",                            & c_eKrasovsky1940,    -43, -163,   45},
	{"Ain El Abd 1970",               "Ain el Abd 1970",                    & c_eInternational,   -150, -251,   -2},
	{"Anna 1 Astro 1965",             "Anna 1 Astro 1965",                  & c_eAustralian,      -491,  -22,  435},
	{"Arc 1950",                      "Arc 1950",                           & c_eClarke1880,      -143,  -90, -294},
	{"Arc 1960",                      "Arc 1960",                           & c_eClarke1880,      -160,  -8,  -300},
	{"Ascension Island 1958",         "Ascension Island 1958",              & c_eInternational,   -207,  107,   52},
	{"Astro B4 Sor.Atoll",            "Astro B4 Sorol Atoll",               & c_eInternational,    114, -116, -333},
	{"Astro Beacon 1945",             "Astro Beacon 1945",                  & c_eInternational,    145,   75, -272},
	{"Astro Pos 71/4",                "Astro DOS 71/4",                     & c_eInternational,   -320,  550, -494},
	{"Astronomic Stn. 1952",          "Astronomic Stn 1952",                & c_eInternational,    124, -234,  -25},
	{"Australian Geodetic 1966",      "Australian Geodetic 1966",           & c_eAustralian,      -133,  -48,  148},
	{"Australian Geodetic 1984",      "Australian Geodetic 1984",           & c_eAustralian,      -134,  -48,  149},
	{"Austrian",                      "Austrian",                           & c_eBessel1841,       594,   84,  471},
	{"Bellevue (IGN)",                "Bellevue (IGN)",                     & c_eInternational,   -127, -769,  472},
	{"Bermuda 1957",                  "Bermuda 1957",                       & c_eClarke1866,       -73,  213,  296},
	{"Bogota Observatory",            "Bogota Observatory",                 & c_eInternational,    307,  304, -318},
	{"Campo Inchauspe",               "Campo Inchauspe",                    & c_eInternational,   -148,  136,   90},
	{"Canton Island 1966",            "Canton Astro 1966",                  & c_eInternational,    298, -304, -375},
	{"Cape",                          "Cape",                               & c_eClarke1880,      -136, -108, -292},
	{"Cape Canaveral mean",           "Cape Canaveral",                     & c_eClarke1866,        -2,  150,  181},
	{"Carthage",                      "Carthage",                           & c_eClarke1880,      -263,    6,  431},
	{"CH-1903",                       "CH-1903",                            & c_eBessel1841,       674,   15,  405},
	{"Chatham 1971",                  "Chatham 1971",                       & c_eInternational,    175,  -38,  113},
	{"Chua Astro",                    "Chua Astro",                         & c_eInternational,   -134,  229,  -29},
	{"Corrego Alegre",                "Corrego Alegre",                     & c_eInternational,   -206,  172,   -6},
	{"Djakarta (Batavia)",            "Djakarta (Batavia)",                 & c_eBessel1841,      -377,  681,  -50},
	{"DOS 1968",                      "DOS 1968",                           & c_eInternational,    230, -199, -752},
	{"European 1950",                 "European 1950",                      & c_eInternational,    -87,  -98, -121},
	{"European 1950 France",          "European 1950 (Mean France)",        & c_eHayford,          -87,  -96, -120},
	{"European 1950 Spain, Portugal", "European 1950 (Spain and Portugal)", & c_eInternational,    -84, -107, -120},
	{"European 1979",                 "European 1979",                      & c_eInternational,    -86,  -98, -119},
	{"Easter Island 1967",            "Easter Island 1967",                 & c_eInternational,    211,  147,  111},
	{"Egypt",                         "Egypt",                              & c_eInternational,   -130, -117, -151},
	{"Finland Hayford",               "Finland Hayford",                    & c_eInternational,    -78, -231,  -97},
	{"France NTF IGN",                NULL,                                 & c_eClarkeIGN,       -168,  -60,  320},
	{"Gandajika Base",                "Gandajika Base",                     & c_eInternational,   -133, -321,   50},
	{"Geodetic Datum 1949 NZ",        "Geodetic Datum 1949",                & c_eInternational,     84,  -22,  209},
	{"Guam 1963",                     "Guam 1963",                          & c_eClarke1866,      -100, -248,  259},
	{"GUX 1 Astro",                   "GUX 1 Astro",                        & c_eInternational,    252, -209, -751},
//rf@	{"Hermannskogel", 682, -203,  480)},
	{"Hermannskogel",                 "Hermannskogel",                      & c_eBessel1841,       653, -212,  449},
	{"Hjorsey 1955",                  "Hjorsey 1955",                       & c_eInternational,    -73,   46,  -86},
	{"Hong Kong 1963",                "Hong Kong 1963",                     & c_eInternational,   -156, -271, -189},
	{"Hu-Tzu-Shan",                   "Hu-Tzu-Shan",                        & c_eInternational,   -634, -549, -201},
	{"Indian Bangladesh",             "Indian Bangladesh",                  & c_eEverestIndia,     282,  726,  254},
	{"Indian India and Nepal",        NULL,                                 & c_eEverestIndia56,   295,  736,  257},
	{"Indian 1954 Thailand",          "Indian Thailand",                    & c_eEverestIndia,     217,  823,  299},
	{"Indian 1960 Vietnam (near 16N)", NULL,                                & c_eEverestIndia,     198,  881,  317},
	{"Israeli",                       "Israeli",                            & c_eClarkePalestine, -235,  -85,  264},
	{"ISTS 073 Astro 1969",           "ISTS 073 Astro 1969",                & c_eInternational,    208, -435, -229},
	{"Ireland 1965",                  "Ireland 1965",                       & c_eAiryMod,          506, -122,  611},
	{"Johnston Island 1961",          "Johnston Island",                    & c_eInternational,    191,  -77, -204},
	{"Kerguelen Island",              "Kerguelen Island",                   & c_eInternational,    145, -187,  103},
	{"L.C. 5 Astro",                  "L.C. 5 Astro",                       & c_eClarke1866,        42,  124,  147},
	{"Liberia 1964",                  "Liberia 1964",                       & c_eClarke1880,       -90,   40,   88},
	{"Luzon Mindanao",                "Luzon Mindanao",                     & c_eClarke1866,      -133,  -79,  -72},
	{"Luzon Philippines",             "Luzon Philippines",                  & c_eClarke1866,      -133,  -77,  -51},
	{"Mahe 1971",                     "Mahe 1971",                          & c_eClarke1880,        41, -220, -134},
	{"Marco Astro",                   "Marco Astro",                        & c_eInternational,   -289, -124,   60},
	{"Massawa",                       "Massawa",                            & c_eBessel1841,       639,  405,   60},
	{"Merchich",                      "Merchich",                           & c_eClarke1880,        31,  146,   47},
	{"Midway Astro 1961",             "Midway Astro 1961",                  & c_eInternational,    912,  -58, 1227},
	{"Minna",                         "Minna",                              & c_eClarke1880,       -92,  -93,  122},
	{"NAD 27 Alaska",                 "NAD27 Alaska",                       & c_eClarke1866,        -5,  135,  172},
	{"NAD 27 Bahamas",                "NAD27 Bahamas",                      & c_eClarke1866,        -4,  154,  178},
	{"NAD 27 Canada",                 "NAD27 Canada",                       & c_eClarke1866,       -10,  158,  187},
	{"NAD 27 Canal Zone",             "NAD27 Canal Zone",                   & c_eClarke1866,         0,  125,  201},
	{"NAD 27 Caribbean",              "NAD27 Caribbean",                    & c_eClarke1866,        -7,  152,  178},
	{"NAD 27 Central",                "NAD27 Central",                      & c_eClarke1866,         0,  125,  194},
	{"NAD 27 CONUS",                  "NAD27 CONUS",                        & c_eClarke1866,        -8,  160,  176},
	{"NAD 27 Cuba",                   "NAD27 Cuba",                         & c_eClarke1866,        -9,  152,  178},
	{"NAD 27 Greenland",              "NAD27 Greenland",                    & c_eClarke1866,        11,  114,  195},
	{"NAD 27 Mexico",                 "NAD27 Mexico",                       & c_eClarke1866,       -12,  130,  190},
	{"NAD 27 San Salvador",           "NAD27 San Salvador",                 & c_eClarke1866,         1,  140,  165},
  	{"NAD 83",                        "NAD83",                              & c_eGRS80,              0,    0,    0},
	{"Nahrwn Masirah Ilnd",           "Nahrwn Masirah Ilnd",                & c_eClarke1880,      -247, -148,  369},
	{"Nahrwn Saudi Arabia",           "Nahrwn Saudi Arabia",                & c_eClarke1880,      -231, -196,  482},
	{"Nahrwn United Arab",            "Nahrwn United Arab",                 & c_eClarke1880,      -249, -156,  381},
	{"Naparima BWI",                  "Naparima BWI",                       & c_eInternational,     -2,  374,  172},
	{"NGO 48",                        "NGO1948",                            & c_eBesselNorway,     315, -217,  528},
	{"Norsk",                         "Norsk",                              & c_eBesselNorway,     278,   93,  474},
	// Algeria
	{"North Sahara 1959",             NULL,                                 & c_eClarke1880,      -186,  -93,  310},
	{"Observatorio 1966",             "Observatorio 1966",                  & c_eInternational,   -425, -169,   81},
	{"Old Egyptian",                  "Old Egyptian",                       & c_eHelmert,         -130,  110,  -13},
	{"Oman",                          "Oman",                               & c_eClarke1880,      -346,   -1,  224},
	{"OSGB 36",                       "Ord Srvy Grt Britn",                 & c_eAiry1830,         375, -111,  431},
	{"Pico De Las Nieves",            "Pico De Las Nieves",                 & c_eInternational,   -307,  -92,  127},
	{"Pitcairn Astro 1967",           "Pitcairn Astro 1967",                & c_eInternational,    185,  165,   42},
	{"Provisional South American 1956 mean", "Prov So Amrican 1956",        & c_eInternational,   -288,  175, -376},
	{"Provisional South Chilean 1963", "Prov So Chilean 1963",              & c_eInternational,     16,  196,   93},
	{"Puerto Rico",                   "Puerto Rico",                        & c_eClarke1866,        11,   72, -101},
	{"Porto Santo 1936",              NULL,                                 & c_eInternational,   -499, -249,  314},
//	{"Potsdam",                       NULL,                                 & c_eBessel1841,       586,   87,  409},
	{"Potsdam Rauenberg DHDN",        "Potsdam Rauenberg DHDN",             & c_eBessel1841,       606,   23,  413},
	// Based on the GOST R 51794-2001 of Russian Federation
	{"Pulkovo 1942 (GOST R 51794-2001)", "Pulkovo 1942 (1)",                & c_eKrasovsky1940, 23.92f, -141.27f, -80.9f, 0, -0.35f, -0.82f, -0.12e-6f},
	{"Pulkovo 1942 Czechoslovakia",   "Pulkovo 1942 CZH",                   & c_eKrasovsky1940,     26, -121,  -78},
	{"Pulkovo 1942 Albania",          "Pulkovo 1942 ALB",                   & c_eKrasovsky1940,     24, -130,  -92},
	{"Pulkovo 1942 Hungary",          "Pulkovo 1942 HUN",                   & c_eKrasovsky1940,     28, -121,  -77},
	{"Pulkovo 1942 Kazakhstan",       "Pulkovo 1942 KAZ",                   & c_eKrasovsky1940,     15, -130,  -84},
	{"Pulkovo 1942 Latvia",           "Pulkovo 1942 LAT",                   & c_eKrasovsky1940,     24, -124,  -82},
	// from http://www.syryjczyk.krakow.pl/Uklad%201965%20i%20inne_T.htm
	{"Pulkovo 1942 Poland",           "Pulkovo 1942 POL",                   & c_eKrasovsky1940, 23.74f, -123.83f, -81.81f},
	{"Pulkovo 1942 Romania",          "Pulkovo 1942 ROM",                   & c_eKrasovsky1940,     28, -121,  -77},
	{"Qornoq",                        "Qornoq",                             & c_eInternational,    164,  138, -189},	
	{"Qatar National",                "Qatar National",                     & c_eInternational,   -128, -283,   22},
	{"Rijksdriehoeksmeting",          "Rijksdriehoeksmeting",               & c_eBessel1841,       593,   26,  478},
	{"Rome 1940",                     "Rome 1940",                          & c_eInternational,   -225,  -65,    9},
	// from http://www.lm.se/geodesi/refsys/eng/refsys-eng.htm
	{"RT90",                          "RT 90",                              & c_eBessel1841,      +414.1055246174f, +41.3265500042f, +603.0582474221f, +0.8551163377f, -2.1413174055f, +7.0227298286f, 0},
	{"RD (Dutch)",                    NULL,                                 & c_eBessel1841,       593,   26,  478},
	{"S-JTSK",                        "S-JTSK",                             & c_eBessel1841,       589,   76,  480},
	{"Santo (DOS)",                   "Santo (DOS)",                        & c_eInternational,    170,   42,   84},
	{"Sapper Hill 1943",              "Sapper Hill 1943",                   & c_eInternational,   -355,   16,   74},
	{"Schwarzeck",                    "Schwarzeck",                         & c_eBesselNamibia,    616,   97, -251},
	{"South American 1969",           "South American 1969",                & c_eSA1969,           -57,    1,  -41},
	{"South Asia",                    "South Asia",                         & c_eFischer1960,        7,  -10,  -26},
	{"Southeast Base",                "Southeast Base",                     & c_eInternational,   -499, -249,  314},
	{"Southwest Base",                "Southwest Base",                     & c_eInternational,   -104,  167,  -38},
	{"Timbalai 1948",                 "Timbalai 1948",                      & c_eEverestTimbalai, -689,  691,  -45},
	{"Tokyo",                         "Tokyo",                              & c_eBessel1841,      -128,  481,  664},
	{"Tristan Astro 1968",            "Tristan Astro 1968",                 & c_eInternational,   -632,  438, -609},
	{"WGS 72",                        "WGS 72",	                            & c_eWGS72,              0,    0,    5},
	{"WGS 84",                        "WGS 84",                             & c_eWGS84,              0,    0,    0},
	{"Viti Levu 1916",                "Viti Levu 1916",                     & c_eClarke1880,        51,  391,  -36},
	{"Wake-Eniwetok 1960",            "Wake-Eniwetok 1960",                 & c_eHough,            101,   52,  -39},
	{"Yacare",                        "Yacare",                             & c_eInternational,   -155,  171,   37},
	{"Zanderij",                      "Zanderij",                           & c_eInternational,   -265,  120, -358},
};

const size_t c_cDatums = sizeof (g_Datums)/sizeof (g_Datums [0]);

CTransformationBase * CDatum::ToWGS84 () const {
	if (pToWGS84 == NULL) {
		if (fWX != 0 || fWY != 0 || fWZ != 0 || fMS != 0)
			pToWGS84 = new CHelmertTransformation    (* pEllipsoid, c_eWGS84, fDX, fDY, fDZ, fWX, fWY, fWZ, fMS);
		else if (fDX != 0 || fDY != 0 || fDZ != 0)
			pToWGS84 = new CMolodenskyTransformation (* pEllipsoid, c_eWGS84, fDX, fDY, fDZ);
	}

	return pToWGS84;
}

bool FindOziDatumToWGS84 (const char * _strOziID, const CTransformationBase * & _pTransformation) {
	if (::strcmp (_strOziID, "Pulkovo 1942 (2)") == 0)
		_strOziID = "Pulkovo 1942 (1)";
	else if (::strcmp (_strOziID, "Pulkovo 1942") == 0)
		_strOziID = "Pulkovo 1942 (1)";

	for (size_t c = 0; c < c_cDatums; ++ c) {
		const CDatum & datum = g_Datums [c];
		if (datum.strOziID && ::strcmp (datum.strOziID, _strOziID) == 0) {
			_pTransformation = datum.ToWGS84 ();
			return true;
		}
	}

	_pTransformation = NULL;
	return false;
}
